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ABSTRACT 


Static and dynamic bifurcations due to nonlinear reactivity interactions in three-temperature 
lumped parameter PHWR dynamical systems are analysed. For the fuel and cool^t temperature 
reactivity feedbacks, we consider both linear and higher-order (quadratic) feedback terms. For 
the moderator, while only linear feedback is considered, we allow both negative and slightly 
positive feedback coefficients. Models based on effective lifetime and prompt-jump 
approximations, as well as those based on one and six groups of delayed neutrons are considered. 
Besides the core and the moderator, two simple models of the steam generator, and 
proportionaland differential controllers are included in the analysis. The control signals may be 
based on fractional change in power and / or coolant temperature. 

Data for analysis is derived from 220 Mwe and 500 Mwe PHWRs as well as from a number of 
other PHWRs in operation outside India. It is found that in all the reactors analysed, Hopf 
bifurcation occurs for negative values of coolant temperature coefficient of reactivity. Since the 
physical values of coolant temperature coefficient of reactivity is positive in PHWRs (under most 
operating conditions), the Hopf bifurcation is unlikely to occur. Even under those operating 

conditions where coolant temperature coefficient of reactivity may be negative, its magnitude is 

-2 

likely to be much smaller than that required for Hopf bifurcation, which is of the order of -10' K' 
* in all the reactors considered. The static bifurcation in these reactors occurs for positive values 
of coolant temperature coefficient of reactivity of the order 10"^ to 10"^ K*^ while the actual 
values is of the order 10'^ to 10"^ K'^ 


ii 



A worst- case analysis based on the most unfavourable combination of dimensionless parameters 
within certain specified ranges is also considered. In this case the values of coolant temperature 
coefficient of reactivity for which a static bifurcation is most likely to occur are found to be of 
the order of +10*^ K’'. While this value is physically realistic for PfTWRs, it must be kept in mind 
that this applies when parameters of all reactors are combined in the worst manner.The worst- 
case parameters for Hopf bifurcation depend strongly on the moderator temperature coefficient of 
reactivity. The critical values of coolant temperature coefficient of reactivity, in particular, range 
from -4 X 10"^ to +2 x 10”^ K"'. 

Inclusion of control reactivity is shown to make bifurcations in the lumped parameter PHWR 
models further unlikely for all negative gains. 
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Chapter 1 


Introduction 

1.1 Motivation 


Nonlinear time evolution and stability of a dynamical system, e.g., an engineering or a physical 
system described by a set of coupled nonlinear ordinary differential equations (ODEs) are 
intimately related to the local bifurcation phenomena that arise in the steady-state form of ODEs. 
Every problem in applications in general contains several physical parameters which may var}’ 
over certain specified sets. Thus it is important to understand the qualitative behavior of the 
system as the parameters vary. A good design for a system will always be such that the 
qualitative behavior does not change when the parameters are varied a small amount about the 
value for which the origmal system was made. However, the behavior may change when the 
system is subjected to large variations in the parameters. A change in the qualitative properties 
could mean a change in the stability of the original system and thus the system must assume a 
state different from the original design. The values of the parameters where this change takes 
place are called bifurcation values. Knowledge of bifurcation values is absolutely necessary for a 
complete understanding of the system. 



The concepts of bifurcation theory, stability theory and nonlinear dynamics are also 
applicable to the dynamics of nuclear systems and can provide cmcial information for the safe 
operation of nuclear reactors. The dynamic models of nuclear systems are nonlinear even with 
linear feedback effects. Nonlinearities in the reactor dynamical systems arise due to strong 
interactions between neutronics, thermal-hydraulics and fission product poisoning. The nonlinear 
model is normally linearised around steady state operating point at full power leadmg to the 
representation of the system by a set of linear ordinary differential equations. These equations are 
transformed into state space form to yield a linear time invariant system for the reactor, which is 
then studied using classical methods applicable to the linear systems. At present time, it appears 
that the techniques of nonlinear systems research have not been as fully and widely utilized in the 
nuclear energy field as in the other branches of scientific knowledge. This work applies some of 
these techniques to lumped parameter PHWR dynamical systems. 

1 .2 Objectives 

The overall objective of the present work is to smdy several fission reactor dynamical systenas, 
representative of lumped parameter PHWRs in a comparative manner, with a view to understand 
the role of intrinsic reactor core nonlinearities arising due to thermal-hydraulic and fission 
product reactivity feedbacks. For this purpose, we choose phenomenological models with 
varying degrees of approximation for the core, moderator, steam generator and the control 
system. 'Bie emphasis will be on the local bifurcation phenomena near steady state operating 
point. Both the static and dynamic bifurcations are studied and the parameter ranges over which 
bifurcations occur are identified. 
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1 .3 Review of Literature 


An extensive review of the literature related to the bifurcation studies applied to PWRs and 
BWRs has been presented by Kalra and Sriram (1998); Manmohan (1996). A brief review of 
these studies applicable to PHWRs is given here. The mathematical and physical importance of 
the nonlinear form of reactor kinetics equations, when the effect of temperature is included, was 
pointed out by Chemick (1951). Robinson (1954, 1955) extended the treatment of Chemick 
(1951) in regard to concept of reactor stability and included the effect of delayed neutrons. The 
concepts of orbital, secular and asymptotic stability are examined in the phase plane 
representation and techniques are pointed out whereby the effect of the delayed neutrons on 
reactor kinetics may be taken into consideration. Hsu (1968), Kastenberg (1968), and Hsu and 
Sha (1969) considered the stability problem for spatially dependent nonlinear reactor systems 
and pointed out that in a practical reactor the delayed neutrons must be included and more than 
one feedback should be considered. Enginol (1985) has analyzed asymptotic stability of nuclear 
reactors with arbitrary nonlinear feedback. Ward and Lee (1987a, 1987b) have used a singular 
perturbation method for the analysis of large amplitude power oscillations in nuclear reactors. 
Adebiyi and Harms(1989, 1990) examine some aspects of linear and nonlinear nuclear reactor 
kinetics by general topological methods. Bergdahl, et al. (1989) analyze the stability and other 
safety-related problems at Forsmark-1 and point out that the instability is caused by the dynamic 
coupling between the neutron kinetics and thermal-hydraulics via void reactivity feedback 
Doming(1989) provides a general discussion of strange atrractors in the context of nuclear 
reactors. Hopf bifurcation iu the context of xenon oscillations is discussed by Rizwan-uddin 
(1989). Yang and cho (1992) present expansion methods for finding nonlinear stability domains 
of nuclear reactor models. 

The data of different PHWRs used in this work is based on the data given in Power Reactors 
(1983) and World Nuclear Industry Handbook (1994). Data of the typical features of standard 
220 MWe and 500 MWe Pressurized Heavy Water Reactors (PHWRs) with natural uranium as 
fuel are given in AERB (1997); Bagchi (1994); NPCIL (1992). A three-temperature model has 
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been developed for PHWRs taking moderator feedback into consideration based on the 
discussions presented by BARC (1987); Rastogi (1976); BHEL. Lewins (1978) presented 
detailed analysis of classical control methods applied to linear reactor dynamical systems. He 
also developed a suitable lumped parameter model for the boiler / steam generator, which has 
been adapted for the present work. 


1 .4 Outline of the Present Work 

Later in this chapter in Sections 1.5 and 1.6 , a brief description of PHWRs and reactivitv' effects 
of temperature is presented. In Chapter 2 , lumped parameter models for core, moderator, steam 
generator and control system based on linear reactivity feedback using delayed neutrons 
represented by one or six groups have been formulated. Aspects of higher-order reactivity 
feedback essential for Static bifurcation analysis are included. Simplications based on effective 
lifetime and prompt-jump approximations are also presented, together with the conditions for 
their validity. All the governing equations are presented in a dimensionless form, and incorporate 
an appropriate translation of the state variables required for later analysis. 

In Chapter 3, the phenomena of Hopf bifurcation and static bifurcation (transcritical) bifurcation 
in the fission reactor dynamical systems are investigated. The numerical scheme for locating the 
bifurcation points in different models is discussed , followed by the calculation of center 
manifold at Hopf and static bifurcations. The effect of system parameters on the critical value 
and the stability parameter is also examined. 

Numerical results on the behavior of the dynamical systems at the close vicinity of bifurcation 
points are presented in Chapter 4. A summary of the work done and conclusions reached, 
together with a few observations on the possible extension of the present work are also presented 
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1 .5 A Brief Description of PHWR 


A Pressurized Heavy Water Reactor (PHWR) is a horizontal pressure tube type reactor using 
natural uranium dioxide fuel with heavy water moderator and high pressure heavy water coolant 
(BARC, 1987; Rastogi, 1976). The fuel and the moderator are contained in a horizontal 
cylindrical shell and tube assembly known as calandria. Cold pressed and sintered pellets of 
uranium dioxide (UO 2 ) in zircaloy tubes form the cylindrical fuel elements. The shell of the 
calandria contains heavy water moderator. 

The moderator and the coolant have their own closed circulating systems. The coolant system, 
called the Primary Heat Transport (PHT) system, is a high pressure and high temperature circuit 
with no bulk boiling. The moderator circuit is at low pressure and low temperature. These two 
are insulated from each other in the reactor core by a gas annulus between the calandna and the 
coolant tubes. The actual data for different PHWRs studied in the present work is given later in 
Chapter 2. A schematic diagram of a PHWR is given in Figure 1 . 

1.6 Temperature Effects On Reactivity 

The reactor operation at power causes changes in the temperature of fuel, coolant and the 
moderator. These changes m temperature cause changes in reactivity which in turn affect the 
power. This change in power again changes the temperatures. Thus a feedback process is 
established. Therefore, variation of reactivity with temperature of the fuel, coolant and moderator 
is an important factor in the reactivity control and the safety of the power reactor. 

The rate of change of reactivity with respect to temperature is defined as the temperature 
coefficient of reactivity, and is expressed as K'^ The effect of the fuel temperature on the 
reactivity is the main contributing factor to power coefficient. The latter with coolant temperature 
coefficient is important for safety assessment during a rapid power transient. The reactivity 
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variation with coolant temperature is important in the assessment of safety aspects associated 
with faults in the coolant system. The usefulness of the moderator temperature coefficient lies in 
its potential applicability to fine reactivity control. The three temperature coefficients of 
reactivity used in the present work are described below. 



Figure 1 : A schematic diagram of a PHWR system with Steam Generator (SG) 
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Fuel Temperature Coefficient 


There are two factors which significantly affect the fuel temperature coefficient: 

(1) Doppler Effect 

When fuel temperature increases, the thermal motion of fuel atoms also increases and the energj' 
range over which resonance capture occurs is broadened resulting in increased resonance 
absorption. This is called Doppler effect. Hence an increase in fuel temperature causes an 

increase in the absorption cross-section of fuel, particularly of • 

(2) Neutron Temperature Effect 

The physical temperature of the fuel affects the fuel neutron temperature because of some 
speeding-up of neutrons in the hot fuel. Due to increase in the temperature of the fuel, the 
average energy of the neutron increases and the spectrum is shifted towards higher energy. 

Coolant Temperature Coefficient 

The effect of coolant temperature on reactivity arises from the changes that occur in the density 
of the coolant and in the neutron speeding-up effect in the fuel as a result of the change in coolant 
temperature, i.e., due to neutron temperature effect. With increase in coolant temperature, the 
density of the coolant decreases, which causes macroscopic absorption crosssections to decrease 
and mean free paths to increase. This in turn affects the reactivity. 


Moderator Temperature Coefficient 

The variation of reactivity with moderator temperature arises from two effects; 
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(1) The density change, associated with the change in moderator temperature, affects the 
moderating properties of heavy water. 

(2) An increase in moderator temperature also causes the thermal neutron spectrum to be shifted 
towards slightly higher energy range at which the cross-sections are different. 
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Chapter 2 


Lumped Parameter PHWR 
Dynamical Systems 

2.1 Neutron Kinetics 

The point-reactor kinetics equations with one group of delayed neutrons can be written as 
(Lewins 1978, Hetrick 1971) : 

— = C , and (2.1) 

dt A 

— = , (2.2) 

dt A 

where 

N = N{t) is the number of neutrons or average neutron density in the reactor, 

C = C(t) is the number of delayed neutron precursors or average delayed neutron 
precursor density, 
t = time(s), 

A = neutron reproduction or generation time (s), 

P = delayed neutron fraction (dimensionless). 



X = decay constant for delayed neutron precursors (s ‘ ), 
p = reactivity or fractional excess multiplication constant (dimensionless). 

2.2 Power Balances 

The three-temperature feedback model is based on three power balances: one for the energy’ 
contained in the fuel elements of the reactor core, second for the energy stored in the coolant 
within the core volume, and another for the fission energy deposited directly in the moderator. 

The power balance for the fuel elements is given by 

c,^= (2.3) 

where 

Ty = Ty(t) is the average fuel temperature (K), 

T +T 

'^c ~ ‘^c(0= ■ ' is the mean coolant temperature (K), 

7) = coolant inlet temperature (K), 

= coolant exit temperature (K), 

P = P{t) is the reactor power (W), 

Cy = heat capacity of fuel elements in the reactor core (JK“’ ), 
and 

Hy = fuel to coolant heat transfer constant (WK'* ). 

The power balance for the coolant \n the reactor core is given by 

(2.4) 
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where 


C, = heat capacity of the coolant in the reactor core (JK ' ), 

OT c = mass flow rate of coolant through the core (kgs ), 

= specific heat of the coolant (J kg'’K‘’). 

The power balance for the moderator ls given by 

= S P-2m,c.,(T„-T„) (2.5) 

Two types of power balance for the moderator have been considered in the present work: 

(1 ) Constant heat removal rate 

In this model the quantity 2m^c„{T„ -T„,,)is a constant and is equal to 5 Hence Eq. 2.5 
becomes 

C,^ = SP-SP.=S(_P-P,) (2.6) 

at 

(2) Heat removal rate proportional to moderator temperature 

In this model Eq. 2.5 represents the power balance for the moderator with the moderator 
temperature T^, appearing on both sides of the equation. 

The symbols in Eqs. 2.5 and 2.6 are as follows : 

= heat capacity of the moderator in the calandria (JK‘^) , 

5 = fraction of fission energy deposited directly in the moderator, 

= flow rate of the moderator through the calandria (kgs'^X 
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c„ = specific heat of the moderator (J kg'^K’*), 

T +T 

'^m ~ ^m(0 = ^ is the mean moderator temperature (K), 

= inlet moderator temperature to the calandria (K), 

= exit moderator temperature from the calandria (K), and 
= steady state power (W). 


The power balance for the boiler / Steam Generator{SG) is given by 




(2.7) 


Two types of power balance for the SG have been derived in the present work : 


(1 ) Constant secondary side temperature 


dT 

In this model the quantity = 0 and hence RHS of the Eq. 2.7 becomes 

= HXT, - rj = 2m, c,{T, - T,) (2.8) 

The second equality in the above equation is based on the assumption that no heat loss takes 
place in the coolant pipes which transport the primary coolant from the reactor to the SG. 


(2) Constant power removal from the secondary side 

In this model the second term on the RHS of Eq. 2.7 is considered constant, i.e., 

mstJi = Po 

In this case the power balance of SG is given by 




(2.9) 
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In Eqs. 2.7 to 2.9, the symbols are as follows : 

Cs — effective heat capacity of the secondary side (JK”' ) 

= steam production rate on the secondary side (kgs 
Ah = specific enthalpy change on the secondary side (Jkg'*), 

= primary to secondary side heat transfer constant (WK'*), 

is saturation temperature corresponding to the secondary side. 

2.3 Steady State Relations 


The steady state operating values of the state variables in Eqs. 2. 1-2.9 will be denoted by 
respectively, and the steady state value of the reactivity is taken as - By 
equating Eqs. 2. 1-2.5 and 2.9 to zero, we obtain 


Pc = 


C = — P 
A 


H 


f rp _p » 

^ Jo ^co 


2m.c.= — 

C C rjp 


( 2 . 10 ) 


^ CO 10 


^ T, , TT 

2m„c„= — —P„ ,and H = 


T -T 

mo mi 


'T 'T 'T' 'T' 

^co -^so ^SO 


where Tjo is the steady state value of T;, the coolant inlet temperature to the reactor. In models 
without the steam generator, T, is considered a constant. 


2.4 Linear Reactivity feedback 


The reactivity p appearing in Eq. 2.1 can, in general, be written as 
P - Po Pjf)'^ Pc Pfi ■ 
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As seen above, if there is no external source of neutrons in the reactor, as is the case considered 
throughout the present work, = 0. For the autonomous system we treat (?) = 0. The control 
reactivity is considered separately in Section 2.10. Further, the dynamical systems developed 
in this chapter are based on linear reactivity feedback which for the three temperature PHWR 
model can be written as 


/? = (r^ - + a, (r, - J (2. 1 1 ) 

where 

c? / = Doppler or fuel temperature coefficient of reactivity (K~' ), 

a ^ = coolant-moderator temperature coefficient of reactivity (K”' ), and 
= moderator temperature coefficient of reactivity (K‘^). 

Reactivity effects of xenon and Higher-order feedback are considered in Section 2.7 and 2.9. 

It is the substitution of Eq. 2.11 in Eq. 2.1 which is one source of nonlinearity in the fission 
reactor dynamical systems. It leads to product terms of P(?),7).(?),r^(?)and TJf) . Instead of 

presenting the resulting equations, we first choose suitable non-dimensional forms of dependent 
state variables and the independent time variable. 

2.5 Dimensionless Equations 


In order to obtain the previous equations in a dimensionless form we choose the following 
translation and scaling of the state variables and time; 


Po N 





Jo 


Tfo 


T/ 


C-C 


T -T 

C*J CO 

•O- — « 

T ~T 

■‘■CO ‘d 


( 2 . 12 ) 
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and 


where 


T -T 
T -T ’ 

CO so 



A 


T -T 

mo 

T -T 

mo mt 


Tj = 7] = Tjo if the SG models are not included in the dynamical systems and 
= if the SG models are included in the dynamical systems. 


The translations of the state variables in the above equations shift the operating or equilibrium 
point to the origin. This is essential for the analysis presented in later chapters. The 
normalizations used could be chosen differently, but the above choice puts the dynamical system 
in the simplest form. Substitution of Eqs. 2.11 and 2.12 into Eqs. 2. 1-2.9 and using Eqs. 2.10 
where needed leads to the following dimensionless system; 


dp 

dt 


= -p + 


ci+aj^+ aX + +a^pTy+a^pTx = 


(2.13) 


d c 
dr 


= b(p-<z). 


dr 


= pO- - &)P -p^f+p^c . 


(2.14) 

(2.15) 


If the secondary side temperature on the SG is constant then we have for the dimensionless 
equation of coolant (model 1, Eq. 2.8) as 


_ pc 

dr d 


( 3 ,- 3 ,), 


(2.16) 


On the other hand, if power removal in the SG is constant (model 2, Eq. 2.9) then the 
dimensionless equation for the coolant is given by 


dr 


^(3^ -3, +(1-^)3,) 


(2.17) 
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For the moderator model 1 (Eq. 2.6), the dimensionless form of the moderator is given by 


^n, 
m 

dr 


= (1P 


(2.18) 


For the moderator model 2 (Eq. 2.8), the dimensionless equation for the moderator is given by 


dr 


9(p-3„,). 


(2.19) 


This completes the dynamical system if steam generator is excluded or model 1 is taken for the 
SG. However, if model 2 of the SG is taken, the above equations are supplemented by the 
following equation governing the dimensionless secondary side temperature: 


dr 




In Eqs. 2.13 to 2.20, the nondimensional parameters are as follows : 


( 2 . 20 ) 


<3, 


P = 




/? 

AP<, 

pCf{T,^-TS 


C, 


r = 




CMTco-Tso) 


^ ^ciT^o-T,) 

P 


<1 = 






( 2 . 21 ) 
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where as mentioned earlier, = T,o if the steam generator is excluded from the dynamical 
system, and Td = Tjo otherwise. 

It can be seen that the original system of the three-temperature linear feedback model has been 
transformed to a dimensionless system with the nine non-dimensional parameters. The 
parameters Qj, and <3,„can be interpreted as the dimensionless measures of temperature 
coefficients of reactivity. The parameter p may imderstood to measure the steady state power 
PgOf the reactor, the other terms in the expression being relatively fixed or showing lesser 
variation for a given reactor. The parameters b, 9, c, g,and rare viewed as simple 
dimensionless combinations of input quantities. 

Considering model 2 of the moderator and the SG (for illustration), the above dimensionless 
equations can be written in the form of a state equation : 

X = V(X,11) 

= U(T1)X + W(X.T1) (2.22) 

where the overdot denotes the derivative with respect to dimensionless time r , X is the state 
vector , p is the set of dimensionless parameters (Eqs 2.21), and V(X,ti) is the vector field whose 
linear part is U('n)X and the nonlinear part is W(X,ti), i.e.. 
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a„,b,p,d,c,q, 



where = 0 is a solution of V(X,ti) = 0 for all ti. It is clear that 
WCO^-n) = 0 and D W(0 ,ti) = 0. 


2.6 Six Groups of Delayed Neutron Precursors 


Six groups appear to be the optimum for an accurate representation of delayed neutron 
precursors. For the calculation of the critical value of a parameter at which a bifurcation may take 
place, it is therefore desirable to use six groups of delayed neutron emitters. This is achieved by 
replacing Eq. 2.14 by six differential equations governing the populations of different delayed 
neutron precursors. 
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i= 1-, 6, 


(2.23) 


where 


C, = population of the precursor group i, 

C,„ = steady state population of the precursor group i. 



A, A 

p ■ 




= delayed neutron fraction for the i precursor group, 
A, = decay constant of the i'* precursor group (s'^), 


/ 


and other variables or constants are as de fi ned earlier in this chapter. Simultaneously, the term 
c in Eqs. 2.13 is replaced by: 






P. . 


where is the relative yield for the precursor group i. For the purpose of parameter 

p 

variations in the six group models, instead of treating each b, as a. variable parameter, we write 



and 




1 



where b has been previously defined (Eqs. 2.21) and is treated as a possibly variable parameter as 
in one group models. 
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2.7 Reactor Models with Xenon Poisoning 


Fission-product poisoning plays a vital role in the operation of power-producing thermal 
reactors. Because of their large thermal-neutron absorption cross sections, two nuclides are of 
particular interest, xenon-135 and samarium-149. Since samarium-149 has a somewhat lower 
cross section for thermal neutrons , it does not contribute as much to the poisoning of a reactor as 
does the xenon 135 and hence it is neglected. 

Denoting by / and X the time-dependent concentrations of I and Xe respectively, we have 

(Lewins, 1978): 

^ = Y I'Lf <j>- ^,1, (2-24) 


dt 




(2.25) 


where 

Y j = yield of iodine atoms per fission, 

Yx~ yield of xenon atoms per fission, 

X, - decay constant for iodine, s"' , 

X^= decay constant for xenon, s , 

<Jx = microscopic absorption cross section of xenon for thermal neutrons, m^ , 

E/ = macroscopic fission cross section, m"’ , and 
(j) = flux of thermal neutrons, m . 

As in the previous sections of this chapter, in order to put Eqs. 2.24 and 2.25 in a dimensionless 
form, we begin by defining 
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and 


(2.26) 


1 = 


I-h 



X = 


x-x„ 




where and X^ axe the steady state concentrations of iodine and xenon obtained by equating 
the right-hand sides of Eqs. 2.24 and 2.25 to zero: 


and 

Aj 


x„ = 




(2.27) 


where T = y j + Y ^ is the steady state value of thermal neutron flux. 


Introducing Eqs. 2.26 into Eqs. 2.24 and 2.25 using (in the context of a point reactor model) 


<Po 


P-Pq 

Po 




we obtain the following dimensionless equations: 


dr 

dr 


= <'(P-I)> and 

= ^(1 - r(Po)P + ^79 - ^(PoX-^{(Po - P)PX > 


where r is the dimensionless time as defined earlier, and 


(2.28) 

(2.29) 




P ’ 



(2.30) 


9^0 + and Y 


-7i 


are the four dimensionless parameters governing the behavior of xenon and iodine. 


The last term in Eq. 2.29 represents the no nlin earity in the governing equation for xenon. In 
order to consider the effect of xenon poisoning on any of the dynamical systems presented 
earlier, the following procedure is adopted: 
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(i) Append Eqs 2.28 and 2.29 to the dynamic model under consideration. 


(ii) Include the reactivity effect due to xenon on the right hand side of the reactivity 
feedback Eqs. 2.11. 

The reactivity feedback due to xenon is given by 


Px=- 


<yxiX-X,) 
vL , 


(2.31) 


where o is the average number of neutrons produced per fission and the other variables are as 
defined earlier. The constants appearing in Eq. 2.31 can be replaced in terms of the 
dimensionless variables defined previously. This yields: 


$x = (2.32) 

where 


= 


r (p„-l 

(Po 


(2.33) 


and X is given by Eq. 2.26 . This reactivity as other feedback reactivities due to (three) 
temperature and effects, interacts parametrically with neutron population when introduced in Eq. 
2.1 yielding an additional linear and nonlinear term in the equation governing neutron 
population. It is straight forward to see that if the resulting equations are collected in the form of 
the state equation (Eq. 2.22), we have the following variants of the dynamical system (with one 
group of delayed neutrons with models 2 of moderator and SG) : 
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PHWR Model with Xenon 
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and 


V= ’'b,P,c,0,q,r,(Po^ 


is the new set of parameters. Here the dimensionless quantities F and y are treated as 

basic physical constants and not as variable parameters. The quantities v and ]3 which appear in 
the expression for are also considered relatively fixed. This leaves as the only new 
nondimensional parameter which may varied due to operation at different steady state neutron 
flux (f)„. 
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2.8 Simpler Dynamical Models 


In the dynamical models presented in the previous sections, the delayed neutrons have been 
treated using one or six precursor groups. One group is convenient and often satisfactory. Further 
simplification is afforded by the effective lifetime model and the prompt-jump approximations. 
In the following subsections these models are briefly described and the governing equations are 
presented in the dimensionless form used in the present work. 


Effective Lifetime Model 

The effective lifetime of neutrons may be defined as (Hetrick, 1971; Lewins, 1978): 

A'=A + ^«^ (2.34) 

XI 

If, in the preceding equations, A is replaced by A^, and the state variable as well as the 
equation(s) corresponding to the delayed neutron precursors are dropped, we obtain the lumped 
parameter dynamic equations based on the effective lifetime of neutrons. This model is 
sometimes easier to analyze but is applicable only when the reactivity and the rate of change 
reactivity are small : 


X 



« \A 


(2.35) 


or, in terms of the dimensionless variables used in this work. 


$ 1 « 1, and 


dr 


«bi$l 


(2.36) 
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where $ — — and b is defined earlier in this chapter. 

H 

Prompt-Jump Approximation 


The prompt-jump approximation or , as it is sometimes called, the zero-lifetime approximation, 
described in detail by Hetrick (1971). Formally, an expansion of the neutron population or 
density in powers of the small parameter A is attempted, the procedure being analogous to the 
method of singular perturbations. If this approximation is used, the point-reactor kinetics Eqs. 2.1 
and 2.2 are replaced by the following single equation : 


dt 


Xp + 


dt 


/3-p 


N 


(2.37) 


In Eq. 2.37 the reactivity p must be substituted through Eq. 2.1 1, in terms of the other variables. 
This makes the Eq. 2.37 for the neutron population nonlinear, as is the case in other models. The 
criterion for the validity of prompt-jump approximation is (Hetrick, 1971): 


1 

dN 

1 

\ 

\ 

N 

dt 

A 


(2.38) 


It should be noted that the criterion (2.38) can only be satisfied if p < , i.e., $ < 1 . In terms of 

the dimensionless variables defined earlier, the prompt-jump approximation can be written as 




(1 + ^)($+ 6 $) 
___ 


(2.39) 


which is applicable if the following condition, equivalent to (2.38), is satisfied: 
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(2.40) 


P 

(l + p)(l-$) 

where the overdot denotes derivatives with respect to the dimensionless time r defined earlier. 

2.9 Higher-order Feedback 

Higher-order feedback terms do not play a significant role in Hopf bifurcation (Manmohan 19981 
However, as will be seen later in this work, they play a crucial role in the study of static 
bifurcations. The higher-order terms of moderator are not considered in the present work. We 
describe below the higher-order feedback reactivity as follows. 


P = Pfb = Pfb“’ + Pfb® + (2-41) 

with 

-r^)+a,(r.-r„) , (2.42) 

+ (2.43) 


where ap and olq denote the second order reactivity feedback coefficients for the fuel and the 
coolant , respectively, and the other symbols are as used earlier in this chapter. It is to be noted 
that the units of ap and are K’^, while those of the linear feedback coefficients, a j and are 

K’^ Introducing Eqs. 2.41-2.43 in Eq 2.1 , and making use of the nondimensional variables and 
parameters already defined, we can write the first equation m reactor kinetics as follows. 

d. ^ 

-^ = -p+(Z^apj 

+ap^j^ (2.44) 
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where 

a,=a,(Tj,-T,yi% 

and all other parameters are as defined earlier. 

2.10 Control Reactivity Feedback 

In this present work we consider a simple control model based on proportional and differential 
controllers. The control may be based on the fractional change in reactor power and / or a 
suitably defined fractional change in reactor coolant temperature. Thus the control reactivity 
can be written as 


Pc + 1, -p + K, (3, + f , -p (2.45) 

where 

Kp = (negative) constant (gain) proportional to the fractional change in power represented 
by the dimensionless power p , 

= (negative) constant (gain) proportional to the fractional change in coolant 
temperature represented by dimensionless 3^ , and 
tp,t^ = time constants associated with the differential controller based on the rate of change 
of power or coolant temperature respectively. 

Inclusion of control reactivity changes only the first equation in the dynamical systems 
presented earlier. If the differential controller is included then the resulting equation is suitably 
recast by transferring the differential terms or substituting their values in terms of state variable 
from other equations. 
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The resulting equation governing the dimensionless power is given below : 






P -{-\ + kp)p+c:+{aj +kj^d)'2y +(ac + ^c +«3„3„ 

+ (2.46) 


where 


1 t KP t B 

^ =— . =— ,and r, = 


P 


A 


A 


Here only the linear terms on the RHS of Eq. 2.46 have been retained. As we will see later, 
analysis based on these linear terms show that bifurcation (static or dynamic) takes place for very 
large (negative) values of Kp or K^., which are unlikely to occur in practice. As a result the 
higher-order terms on the RHS are not required in the present work. 


2.11 Input Data 


The set of input values that we need to know for each reactor model to carry out the analysis and 
computations have been arrived at in the preceding sections of this chapter. The input data are 
evaluated from the reactor data available from Power Reactors (1983) and World Nuclear 
Industry Handbook (1994). The fuel bundle and other data required for the standard 220 MWe 
and 500 MWe PEIWR are taken from the AERB (1997); Bagchi (1994). Data of about 60 
different PHWRs were grouped into 14 categories based on similarity of their operating and 
design features. Under different operating conditions, the input data can be widely different. The 
data for normal full power operation is summarized in Table 1. 
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Table 1 : Input data* for PHWRs 



To 

Cf 

Cc 

Tfo 

T 

^co 

T 
^ 10 

Reactor 

MWt 

MJK'^ 

MJK"' 

K 

K 

K 

220MWe-India 

801 

17.9 

16.5 

1313.5 


522.0 

SOOMWe-India 

1725 

35.8 

41.3 

1253.5 

555.0 

533.0 

Bruce 1-4 

2606 

38.1 

62.3 

1438.0 

547.5 

522.0 

Bruce 5-8 

2832 

38.1 

67.7 

1438.0 

554.0 

530.0 

Ceraavada 1-5 

2180 

26.7 

52.1 

1386.0 


539.0 

Darlington 1-4 

2774 

38.1 

66.3 

1423.0 

562.5 

538.0 

Embaise-Cordoba 

2103 

27.4 

50.3 

1386.0 

562.0 

539.0 

Gentilly 2 

685 

28.0 

52.1 

1386.0 


539.0 

Kanupp 

433 

9.72 

8.91 

1344.5 

542.5 

519.0 

Pickering 1-8 

1754 

29.7 

41.9 


544.4 

522.0 

Point Lepreau 1 

2156 

28.0 

51.5 

1386.0 

561.5 

■9IM 

■■■■■ 

Wolsong 1-2 

2180 

27.5 

52.1 

1386.0 

561.3 

539.6 

Atucha 1 

1179 

12.4 

24.3 

1559.5 


535.0 

Atucha 2 

2160 

27.2 

51.6 

1559.5 

567.9 



t The meaning of each symbol is given in the Nomenclature and explained in the text. 


Sample calculations of and Q for the 220 MWe are shown in Appendix A. The input data 
for the lumped parameter models of the moderator and steam generators is given in Table 2. A 
sample calculation for Cm is given in Appendix A. The value of is taken to be One Full Power 
Second (FPS). The values of moderator temperatures are as given in Power Reactors (1983) and 
World Nuclear Industry Handbook (1994). The fraction of fission energy directly deposited in 
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the moderator (5) is taken to be 5% in this work. The steady state SG secondary side 
temperatures based on the available data on the secondary side pressures (NPCIL 1992). 


Table 2: Input data for the lumped parameter model of moderator and steam generator 


Symbol 

220MWe-India 

500MWe-India 

Cm 

330 MJK'‘ 

840 MJK'‘ 

5 

0.05 

0.05 

Tmo i 

328 K 

343 K 

Tm. 

318K 

328 K 


-5 X 10'^ K*' 

-5 X 10’" K'^ 

Cs 

1 FPS 

1 FPS’ 

Tso 

524 K 

529 K 

i 


I Xhe meaning of each symbol is given in the Nomenclature and explained in the text. 

Table 3 gives input data, and the corresponding dimensionless parameter values used forone- 
group reactor kinetics. The dimensionless parameters are for the 220 MWe reactor.The coolant 
temperature coefficient a^. ( or the corresponding dimensionless parameter a^) is treated as the 
variable parameter in the PHWR models. ^ 

Tables 4 and 5 give six-group of delayed neutron data and the input constants and dimensionless 
parameters related to xenon and iodine (Lewins, 1978; Hetrick, 1971). The six-group data used is 
for thermal fission of uranium-235. The dimensionless parameters shown in Table 5 are for (t)o = 
10^® m'^s'* and wiU change when a different steady state operating flux is considered. Fianlly 
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Table 6 lists the dimensionless thermal-hydraulic parameters (p,c,0,q,r — all defined earlier) for 
the different categories / groups/ types of PHWRs considered in this work. 

Table 3 : Input data and dimensionless parameters for one-group reactor kinetics 


Input Values 

Dimensionless Parameters 

Symbol 

Value with Units 

Symbol 

Value 

P 

0.006 

b 

0.0075 

A, 

0.075 s‘^ 


-2.64 

A 

6.0 X 10"^ s 


0.022 

ttf 

-2 X 10'^ K‘‘ 


-2.65 

etc 

6xlO^K'‘ 



am 

-5 X 10'^ K‘^ 




Table 4: Input data for six groups of delayed neutrons 


Group 

i 

Relative Yield 

ai 

Decay Constant 
Xi(s-') 

1 

0.038 

0.0127 

2 

0.213 

0.0317 

3 

0.188 

0.115 

4 

0.407 

0.311 

5 

0.128 

1.40 

1 ^ 

0.026 

3.87 
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Table 5 ; Input data and dimensionless parameters for xenon and iodine 


Input Data 

Dimensionless Parameters 
(At V = 2.5, <j)o = 10** m'V*) 




Value 

Y. 

0.064 

i 

2.1742 X 10'^ 


2.87 X 10'^ s'* 

X 

1.5833x 10'^ 

Yx 

0.003 

9 

13.9187 


2.09 X 10'^ s'* 

Y 

0.9552 

<7x 

2.70 X 10'**^ m** 







Table 6 : Dimensionless parameters' for PHWRs based on data given in Table 1 


Reactor 

P 

C 

e 

q 

r 

220MWe-India 

0.005809 

1.08606 

0.027795 

0.001192 

0.005000 

500MWe-India 

0.006891 

0.86884 

0.030534 

0.000685 

0.003846 

Bruce 1-4 

0.007685 

0.61113 

0.027838 

0.000947 

0.005405 

Bruce 5-8 

0.008413 

0.56248 

0.026432 

0.001535 

0.003333 

Cemavada 1-5 

0.009830 

0.51613 

0.025974 

0.000986 

0.003125 

Darlington 1-4 

0.008466 

0.57409 

0.027684 

0.001205 

0.002985 

Embaise-Cordoba 

0.009304 

0.54565 

0.027155 

0.000710 

0.003030 

Gentilly 2 

0.002965 

0.53763 

0.025974 

0.000318 

0.000382 

Kanupp 

0.005555 

1.09090 

0.028468 

0.000551 

0.005405 

Pickering 1-8 

0.006660 

0.70913 

0.024656 

0.000649 

0.006494 

Point Lepreau 1 

0.009329 

0.54395 

0.025414 

0.001012 

0.003077 

Wolsong 1-2 

0.009605 

0.52841 

0.025638 

0.000986 

0.003096 

Atucha 1 

0.009430 

0.51175 

0.016593 

0.001066 

0.004348 

Atucha 2 

0.008000 

0.52730 

0.017245 

0.001724 

0.002571 


t The meaning of each dimensionless parameter is given in the Nomenclature and explained 
in the text. 
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Chapter 3 


Analytical and Computational 
Techniques 

3.1 Operating Point Stability and 
Bifurcation Points 


The steady stateof a linear system is determined by equating the rates of changes X, to zero. 
Inevitably in a real system there will be small disturbances and we study the resulting transients. 
If these disturbances, whatever their inital form, die out, the system is stable. If they grow,or 
even if only some of them grow, the linearity implies that they will go on growing without bound 
\ the system is unstable. There is a possiblity of a behavior on the limit of stability, where an 
oscillatory behavior is maintained whose amplitude neither increases nor decreases from the 
initial disturbance. 

The solution structure of the linearised system 


X = U{ri)X 



associated with the nonlinear dynamical system (Eq. 2.22) is determined by the eigenvalues and 
the corresponding eigenspaces of the Jacobian matrix U(t|). As noted earlier 0 is an 
equilibrium or fixed point for the dynamical system for all values of q. For those values of q for 
which none of the eigenvalues of U(q) is on the imaginary axis, the fixed point is called 
hyperbolic. For these parameter values, the stability properties of the nonlinear system (Eq. 2.22) 
locally coincide with those of the linear system in some neighborhood of q. If for some q= q ’ , 
some or more of the eigenvalues are on the imaginary axis, the equilibrium point is called non- 
hyperbolic. In the vicinity of q ’ the stability properties of the system (Eq. 2.22) can only be 
understood by including the high-order nonlinear terms in the analysis. In such parameter- 
dependent systems, it is important to understand the qualitative behavior of the system as the 
parameter change. A bifurcation describes a qualitative change in the dynamics when the 
parameters in the system are varied. The values of the parameters where this change takes place 
are called critical values or bifurcation values. Knowledge of bifurcation values is absolutely 
essential for the complete understanding of the system. The best understood bifurcations are 
those for which a local analysis is possible. Among these, the most important bifurcations are of 
codimension one, i.e., bifurcations which are essentially one-parameter phenomena and can 
occur m a persistent way in a one dimensional parameter space (Bruno, 1989). In this chapter we 
examine the phenomena of Hopf bifurcation and static bifurcation in the dynamic models of the 
fission reactors of PHWR type, using a local analysis based on the theory of center manifold and 
normal forms (Wiggins, 1989; Guckenheimer and Holmes, 1983). 

Hopf Bifurcation 

A Hopf bifurcation occurs if for some q= q * (the critical or bifurcation value), there exits a pure 
imaginary pair of eigenvalues which crosses the imaginary axis transversally, as q is varied, 
while all the other eigenvalues remain in the left-half plane(Figure 2b). If the minimum 
dimension of the control parameter to be varied is one, we have a Hopf bifurcation of 
codimension one. Bifurcations of higher dimension are of less practical interest. At Hopf- 
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bifurcation, an equilibrium point bifurcates into limit cycles. A hopf bifurcation is also 
sometimes refereed to as a complex or dynamic bifurcation. 

Static Bifurcation 

A static bifurcation occurs if for some t|= ti * a zero eigenvalue crosses the imaginary axis 
through the origin as r\ is varied, while all other eigen values remain in left-half plane(Figure 
2a).Saddle-node, transcritical and pitchfork are three types of static bifurcations. The type of 
bifurcation where on one side of a parameter value there are no fixed points and on the other side 
there are two fixed points is refered to as saddle-node bifurcation. If at the bifurcation point 
exchange of stability occurs then it refers to transcritical bifurcation. In pitchfork bifurcation, the 
operating fixed point which is stable becomes unstable creating two new stable fixed points at 
the bifuracating point. A static bifurcation is referred to as a real bifurcation. 

A Hopf bifurcation is illustrated in Figures 4a,4b and 4c, while the the three types of static 
bifurcations are illustrated in Figure 3a, 3b and 3c. 
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Figure 2 ; A schematic diagram of (a) static bifurcation (b) hopf bifurcation 
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Figure 3: Schematic 


,ic diagrams of (a) saddlo-oode (b) pitchforlc and (c) transcridcal bifincadons 
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Figure 4 : (a) A schematic diagram of a normal (supercritical)Hopf bifurcation 



Figure 4 : (b) A schematic diagram of a subcritical (inverse) Hopf bifurcation 
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Figure 4 : (c) Coexistence of unstable and stable solutions in a subcritical Hopf bifurcation 
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3.2 Normal Forms at Bifurcation Points 


Several analytical and computational techniques are needed to locate the bifurcation points and 
to investigate the behavior of the dynanaical system in the vicinity of these points. These has 
been discussed in detail by Wiggins (1989) , Parker and Chua (1989), . The techniques required 
for the analysis of Hopf Bifurcation are also developed in Manmohan (1996). Below we 
summarize the methods used in the present work for the study of both the Static and Hopf 
Bifurcation in PHWR systems 


Poincare Normal Form 

A systematic way to ascertain the effect of the nonlinear terms is to reduce the dynamical system 
to a two dimensional center manifold, which is then transformed to the Poincare normal form. 
The proof that this can be done may be foxmd in Wiggins (1989). Here, we state the main results 
related to the Poincare normal form and introduce some notations used in the calculations. The 
Poincare normal form in the cartesian coordinates, Xj and X 2 for Hopf bifurcation, can be written 
as 


Xi = a(s)xi - (oi£)x2 + [a{_£)x^ - h{s)xT}{x^^ + x/)+ (3.1) 

X 2 = co{s)x^ + a{s)x 2 +[K^)Xi +a{s)x 2 ]{xi^ +^ 2 ^)+ (3-2) 

where the normal form parameters a{£) and b(s) are defined for all sufficiently small 8. The 
parameter a(0) , called the stability parameter, determines the stability of bifurcating penodic 

solutions (Wiggins 1989): 

( i ) The case a(0) < 0 : The bifurcation is normal (supercritical). 
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( ii ) The case a(0) > 0 : The bifurcation is inverse (subcritical). 


Remark : A simple proof of the above statements, facilitated by considering Poincare normal 
form in polar coordinates, may be found in Wiggins (1989). 

The Poincare normal forms for static bifurcation in Cartesian coordinates are as given below: 

x = £±x^ (saddle-node), 

X = £c + (pitchfork), 

X = sx + x^ (transcritical). 


In the above equations s is the variable parameter. Let p € 77 , and p* be the value of p at which 
the bifurcation occurs. Then s is defined as 


s = 


♦ 


3.3 Location of Bifurcation Points 


Several techniques have been discussed in the literature for the location of the bifurcation point 
(Kubicek and Marek, 1983). In general the algorithms used are iterative in nature and make use 
of the Newton-Raphson method to find the critical value of the parameter. 

The critical value p* of p e 77 is located by solving the equation 

a(p)=0 (3.3) 

where a (p) denotes the real part of a complex conjugate pair of eigenvalues of the matrix U(ti), 
evaluated at the steady state / equilibrium (operating) point. 
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The QR Method 


The eigenvalues of the matrix U(t]) for a given set of parameter values are calculated using the 
QR algorithm. The Jacobian matrices associated with the dynamical systems studied in the 
present work are, in general, nonsymmetric. Since a nonsymmetric matrix may not be balanced, 
it is first replaced by a balanced matrix with identical eigenvalues (Press, et al, 1993). This 
reduces the sensitivity of the eigenvalues to rounding errors during numerical computations. The 
matrix is then reduced to a simpler Hessenberg form and the eigenvalues are calculated by using 
the standard QR algorithm applicable to real Hessenberg matrices. 

Newton-Raphson Iteration 

The Newton-Raphson algorithm is used to find the solution of Eq 3.3. Since the success of a 
locally convergent scheme depends critically on having a good first guess for the solution, it is 
desirable to use a globally convergent method which guarantees some progress towards the 
solution at each iteration. Perhaps the best strategy is to bracket the root and then refine it by a 
combination of the Newton-Raphson and the bisection methods. This keeps the root within the 
brackets while taking advantage of the rapid local convergence of the Newton-Raphson method. 
We have mostly used this kind of approach in this work. 


3.4 Calculation of Center Manifolds 


The center manifold theory provides a means for systematically reducing the dimension of the 
state space which needs to be considered while analyzing the bifurcation of a given type. The 
dynamical system becomes further simplified if we put the center mamfold in a normal form. 
This is based on the idea of introducing successive coordinate transformations to simplify the 
analytic expression of the vector field on the center mamfold. In fact mathematicians have 
proved that under fairly general conditions, the local theory of codimension one bifurcations 
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from a fixed or equilibrium point can be reduced to few archetypes (Ruelle, 1989; Berge, et al, 
1984). It then suffices to decide which of them fits a given problem in order to obtain a 
qualitative picture of the flow in the neighborhood of the bifurcation point. 

The analysis based on the center manifold and normal form cadculations is facilitated if the 
operating or steady state point is translated to the origin. 


Hopf Bifurcation 


Calculation of center manifolds and normal forms requires considerable amount of algebra, 
[referred to as Horrendous by Wiggins (1989) ] particularly for the 2-dimensional Hopf 
bifurcation. These algebraic aspects axe investigated in detail by Manmohan (1996) and in the 
present work we used the same computer program to determine the stability parameter, a(0), 
occurring in the Poincare normal form. We first consider the parameter independent case and let 
U=U(r|*) and W(X) = W(X, q*) to simplify the notation. 

In order to calculate the center manifold of the dynamical system 

X = UX+W(X) (3.4) 

letX= SY 

where S is the matrix whose columns are the eigenvectors of U. 

Y=(S-‘US)Y + S-‘W(SY) 

= AY + S'^W(SY) 

Equation 3.2 can be separated into two parts : 

u = Aw + F(u,v) 
v = Bv + G(u,v) 


Thus we have 

(3.5) 

(3.6) 


(3.7) 


where 
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A r ® *^o~\ fu 

A= . , Y= 

_“®o 0 J [y 

Then tlie introduction of another nonlinear transformation 
V = h('w) 


followed by elimination of quadratic terms in the resulting equations is required to obtain the 
Poincare normal form. For details reference should be made to Wiggins(1989). 
Mamnohan( 1 996). 

Static Bifurcation 

The transformations required for the calculation of center manifold/normeil form for static 
bifurcation are similar to those for the Hopf bifiircation, but considerably simpler due to one- 
dimensional nature of the manifold. We illustrate it using one-group 1lTHU^[paodel., The center 
manifold u, (3 .7) is one-dimensional and we can write it as 

u=F(u,v) 

let Vi = hi(u), i=l,2, n-1. 

1 2 
= biU 

As shown in Wiggins(1989), each hi(u) must satisfy 

^F(u,h(u))-Aihi -Gi(u,h(u)) = 0, i= l,2,...n-l. (3.9) 


where 
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A=. 


0 . ... 0 

•A. ... . 

Aj . . 




is a diagonal matrix. 


It is straightforward to see that if only linear feedback is used solution of Eq. 3.9 yields, 

b, = 0, i=l,...,n-l. 

Thus, the center manifold / normal form has no quadratic terms. If all b,’s are zero, and the 
original nonlinear terms in W(X) do not contain cubic terms (as is the case with linear feedback), 
the center manifold or normal form also do not contain any cubic terms. 


Thus with linear feedback, all the three types of static bifurcation are ruled out, and the behavior 
of the nonlinear system remains divergent as would have been the case in the linear system. This 
is in conformity with the fact that, with linear feedback, no additional equilibrium (fixed) point 
originates from the operating point, i.e., (0,0, ...0). 


If higher-order feedback effect are present, then, as shown in sec-2.9, these introduce additional 
quadratic and cubic terms in the system (Eq. 2.44) . Due to this cubic terms, an additional 
equlibrium solution/fixed point of interest emerges. It is easy to see that without the SG, this 
fixed point is given by (using one group of delayed neutrons) 






and with the SG, the new fixed point is given by 


(3.10) 
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p=c=3„=0, 3 ^=& 3 ^, 3 


ay&+a^+a^ 
Opd^ +ar+a 


> 3,=5c 


(3.11) 


M 


In the nominal calculations presented in Chapter 4, we have taken a„ = 0 
It is straight to see that the bifurcation point corresponds to 


=0 

or a^r^ + a,. = 0 

as the case may be. Thus the new fixed / equilibrium point represented by Eq. 3.14 or 3.15 
originates firom the origin, i.e., from the operating point. Since only one new fixed poini 
bifurcates firom the origin, it indicates transcritical bifurcation, and excluding the possibility of 
saddlenode and pitchfork bifurcations. (The saddle node bifurcation is also excluded firom the 
fact that the original equilibrium point, (0,0,. ..0) does not remain stable beyond the bifurcation 
point). The fact that at the bifurcation point, exchange of stabihty occurs, is brought out in 
diagrams shown in Chapter 4. 

As can been seen firom Eq. 2.44, the higher-order feedback effects introduce into the system not 
only the cubic terms, but additional quadratic terms. As a result, the quadratic terms in the 
normal form, can no longer be zero and thus, the normal form will correspond to that of 
transcritical bifurcation. 


3.5 Calculation of Trajectories 

Calculation of trajectories by numerical integration is perhaps the most important and expensive 
task in the simulation of a dynamical system. It is also a prerequisite for the application of 
geometric tools or mathematical constructs such as Poincare sections, phase portraits and tiie 
calculation of characteristic multipliers or Lyapunov exponents. 
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Over the years, a large number of algorithms have been developed to approximate the behavior 
of a continuos-time system on a digital computer. A proper choice of the method and a careful 
application play a crucial role in the computational accuracy and efficiency. Among the available 
methods, the major practical choices are the single step methods, the multistep predictor- 
corrector methods, and the Bulirsch-Stoer extrapolation methods. The most popular among the 
single step methods are the Runge-Kutta methods. Both explicit and implicit implementations of 
the methods mentioned above are available in literature. The choice depends on the stifftiess of 
the system to be integrated. In the present work an implicit implementation of Bulirsch-Stoer 
method given by Press, etal (1993) is used. 

3.6 Poincare sections and Characteristic 
Multipliers 

In order to simplify the phase space diagrams of a dynamical system of order n, one introduces 
an (n-1) dimensional surface of sections in the phase space and, instead of studying a complete 
trajectory, one monitors only the points of its intersection with this surface. A set of points of 
intersection of the trajectory with a hypersurface, as the trajectory crosses the surface from one 
side to the other, is refeued to as Poincare section. Mapping an intersection point onto the 
subsequent intersection point is referred to as the Poincare (return) map. 

The stability of periodic solutions or closed orbits is determined by the characteristic (Floquet) 
multipliers or roots. These can be regarded as a generalization of the eigenvalues at an 
equilibrium point. According to standard Floquet theory, the stability of a periodic orbit is 
determined by the eigenvalues of the monodromy matrix <I>(T) ( also referred to as Floquet 
transition matrix), which is simply defined as the value of the fundamental matrix at the time 
period T, with the initial value as the umt matrix I (Kubicek and Marek 1983). If these 
eigenvalues lie within the unit circle in the complex plane, the periodic solution is stable. It is 
well known that one of the eigenvalues of d>(T) is always umty (Husseyin, 1986, Parker and 
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Chua 1989). Thus a periodic orbit can never be asymptotically stable. The remaining eigen 
values of 0(T) determine the orbital stability and are called the characteristic multipliers. 

In the present work, we confirm the existence of limit cycles (whenever these occur) using 
Poincare sections, and calculate the characteristic multipliers using a method based on Poincare 
return maps (Manmohan 1996). The numerical results for both the transcritical and Hopf 
bifurcations are reported and discussed in the next chapter. 
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Chapter 4 


Numerical Results, Discussion, and 
Conclusions 


Numerical simulations are required to conform the predictions of the theoretical considerations 
given in Chapter 3. Furthermore, the theory is valid only in an arbitrarily small neighborhood of 
the bifurcation points. As one moves away from this neighborhood, numerical calculations are 
often the only means of obtaining information about the nature of the solutions. In this chapter, 
we present the numerical results for the lumped parameter PHWR dynamical systems developed 
in Chapter 2. 

4.1 Stability at Full Power Operation 


Before proceeding to calculate the critical values of the parameters for which the bifurcations 
may occur and simulating the behavior subsequent to bifurcation , it is desirable to confirm the 
stability of the systems considered at their normal operating values. In reactor dynamical 
systems, this indicated by the reactor period, Tp , or the e-folding time Glasstone (1994), 


( \ 



\dt) 


where cd^ is the real part of the eigenvalue closest to the imaginary axis. The eigenvalues are 
calculated for the Jacobian matrix U, (Eq. 2.22), evaluated at the normal operating point. A 
negative reactor period indicates stability at the point under consideration. 

Table 7 shows the (negative) reactor periods for the reference PHWR (220 MWe) using five 
dynamical models presented in Chapter 2. It can be seen that for all the models ( reactor periods 
are shown for the core, core along with moderator , core with moderator and SG), reactor 
stability is indicated at full power steady state operation. This is also the case for the other 
reactors studied. This is indicated in Table 8 using six groups of delayed neutrons. 

Table 7: Reactor Period for full operation at full power for the Reference PHWR(220 MWe) 


Reactor Model 

Negative Reactor Period (s) 

Core 

Core + Moderator 

Core + Moderator + SG 

No delayed neutrons 

32.7 

81.95 

1.18x10^ 

Effective- lifetime model 

35.4 

81.75 

1.18x10^ 

Prompt-jump approximation 

9.80 

81.73 

8.95 X 10^ 

One group of delayed neutrons 

9.77 

81.73 

1.18x10^ 

Six groups of delayed neutrons 

83.3 

82.53 

1.18x10^ 





'i % 
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Table 8 : Reactor Periods at full power operation for the PHWRs shown in Table 1 
( using six groups of delayed neutrons ) 



Negative Reactor Period* (s) 

Reactor 

Core 

Core + Moderator 

Core + Moderator + SG 

220 MWe-India 

8.33 X 10^ 

8.25 X lO' 

1.18x10" 

500 MWe-India 

8.33x10' 

1.39x10" 

1.11 xlO" 

Bruce 1-4 

8.33x10' 

1.02x10" 

1.38x10" 

Bruce 5-8 

8.33 X lO' 

8.32x10' 

1.37x10" 

Cemavada 1-5 

8.34x10' 

9.71 X lO' 

1.29x10" 

Darlington 1-4 

8.33 X lO' 

8.16x10' 

1.35 X lO" 

Embaise-Cordoba 

8.34x10' 

1.19x10" 

1.29x10" 

Gentilly 2 

8.30x10' 

3.02x10" 

1.37x10" 

Kanupp 

8.33 X lO' 

1.18x10" 

1.22 X lO" 

Pickering 1-8 

8.32x10' 

1.48 X lO" 

1.37 X lO" 

Point Lepreau 1 

8.34 X lO' 

9.45 X lO' 

1.29 X lO" 

Wolsong 1-2 

8.34x10' 

9.71 X lO' 

1.29 X lO" 

Atucha 1 

8.32x10' 

9.04x10' 

1.55 X lO" 

Atucha 2 

8.32x10' 

8.31 X lO' 

1.54x10" 


t A negative reactor period indicates stability at normal operating point of the reactor based on 
the data and reduced -order models used for the calculations. 


The effect of different models of moderator and SG on reactor period is shown in Table 9 and 
unless otherwise specified, the models 2 of moderator and SG (Eqs. 2.5 and 2.9) are used along 
with the core in all the analysis. 
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Table 9 ; Effect of different models of moderator and steam generator (SG) on reactor period for 
the reference reactor (220 MWe) ( normal operation at full power using one group of delayed 
neutrons) 


Moderator Model 

SG model 

Negative Reactor Period (s) 


Constant secondary side 

temperature 

2.64x10^ 

Constant heat removal rate 




Constant power removal in 

theSG 

7.91 X 10^ 


Constant secondary side 

8.17X10^ 

Heat removal rate 

temperature 


proportional to moderator 



temperature 

Constant power removal in 

1.18X10^ 


theSG 



4.2 Critical Parameter Values 


For the purpose of bifurcation studies, we have been treated the coolant temperature coefficient 
of reactivity , , as the variable parameter, keeping all other parameters fixed. ( A simultaneous 

variation of all parameters with a view to find the worst combination is presented later.) The 
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critical values, a „ for the static and dynamic bifurcations for the reference reactor 220 MWe 
core are presented in Table 10 using 


• No delayed neutrons 

• Effective life time model 

• Prompt-jump approximation 

• One group of delayed neutrons 

• Six groups of delayed neutrons 


Table 10 : Critical Values of the coolant temperature coefficient of reactivity ) for static and 
dynamic bifurcation for the Reference PHWR core 


Reactor Model 

Critical Values a g K‘^ 

Static 

Dynamic 

No delayed neutrons 

7.20x10^ 

-8.19x10"* 

Effective- lifetime model 

7.20x10'^ 

-8.48 X 10'^ 

Prompt-jump approximation 

7.20x10'^ 

— 

One group of delayed neutrons 

7.20 X lO"^ 

-6.270 X 10’^ 

Six groups of delayed neutrons 

7.20x10"^ 

-5.428 X 10'^ 


It can be seen that the critical values for Hopf bifurcation considerably depend on the model 
employed. In fact the prompt-jump approximation does not give any critical value for the Hopf 
bifurcation. For this reason, unless otherwise indicated, we have used one-group and six-group 
models only. 
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The effect of moderator and SG on the critical value, a*^ , is brought in Tables 1 1 and 12, for the 
static and Hopf bifurcations respectively, using one and six groups of delayed neutrons. While 
the effect of the number of groups of delayed neutrons is relatively small, the effect of moderator 
and SG on the critical value is considerable. 


Table 1 1 : Effect of moderator and steam generator (SG) on the critical values of a,; for the 
static bifurcation in PHWR (220 MWe) 


Core Model 

Critical Values (static), K'‘ 

. — 

Core 

Core + Moderate] 

Core + Moderator + SG ' 

One group of delayed neutrons 

Six group of delayed neutrons 

7.20 xlO"* 
7.20x10'^ 1 

7.42 X lO"^ 
7.42x10"^ 

2.0 X 10’^ 

2.0 X 10'^ 


It is observed that xenon destroys the possibility of static bifurcation by introducmg a complex 
eigenvalue at the bifurcation point. Hence its effect only on dynamic bifurcation is studied usmg 
PHWR (220 MWe) along with moderator and SG using one and six groups of delayed neutrons. 
This is shown in Table 13. For comparison the critical values without the effect of xenon on 
dynamic bifurcation are also presented in Table 13. It can be seen that the effect of xenon on the 

critical values is negligible. 
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Table 12 ; Effect of moderator and steam generator (SG) on the critical values of for the 
dynamic bifurcation in PHWR (220 MWe) 


Core Model 

Critical Values (dynamic) , K'^ 

Core 

Core + Moderator 

Core + Moderator + SG 

One group of delayed neutrons 

Six group of delayed neutrons 

-6.269 X 10'^ 

-5.442 X 10'^ ' 

-6.2720 X 10’^ 

-5.4303 X 10’^ 

-6.9662 X 10'^ 
-6.1063 X 10'^ 


Table 13: Effect of xenon on the critical values of for dynamic bifurcation in PHWR (220 
MWe) (Core + Moderator + SG) 



Critical Values a VK'^ 

Core Model 

Without Xenon 

With Xenon 

One group of delayed neutrons 
Six group of delayed neutrons 

-6.9662 X 10‘^ 
-6.1063 X 10'^ i 

-6.9640 X 10'^ 

-6.1040 X 10'^ 
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Table 14 shows the critical values for the other reactors using one group of delayed neutrons and 
mcluding the moderator and SG in the model. It can be seen from Tables 10 to 14, that the static 
(transcntical) bifurcation occurs for the positive values of the coolant temperature coefficient of 
reactivity (aj.), while the dynamic (Hopf) bifurcation occurs for the negative values of the same. 
Since the physical value of a^. is positive in PHWRs (under most operating conditions), the Hopf 
bifurcation is unlikely to occur. Even under those operating conditions where may be 
negative, its magnitude is likely to be much smaller than that required for the Hopf bifurcation, 
which is of the order of -10 ^ K ' in all the reactors considered. The static bifurcation in these 
reactors occurs for the positive values of of the order 10'^ to 10"^ K”\ while the actual values is 
of the order 10'^ to 10"^ K"*. 

Table 14: Critical Values of coolant temperature coefficient of reactivity a^for static and 
dynamic bifurcations for the PHWRs listed in Table l(Core +Moderator) using one-group of 
delayed neutrons 



Critical Values a*c K'^ 

Reactor 

Static 

Dynamic 

220MWe-India 

7.22 

X 

10"^ 

-6.2720 X 10’^ 

500MWe-India 

6.86 

X 

lO"^ 

-5.2701 X 10'^ 

Bruce 1-4 

7.48 

X 

10"^ 

-4.1681 X 10’^ 

Bruce 5-8 

7.78 

X 

10"^ 

-4.0991 X 10'^ 

Cemavada 1-5 

8.04 

X 

10"^ 

-3.9009 X 10'^ 

Darlington 1-4 

7.53 

X 

10"^ 

-3.9867 X 10'^ 

Embaise-Cordoba 

7.69 

X 

10"^ 

-3.9072 X 10"^ 

Gentilly 2 

7.89 

X 

lO"* 

-1.0634x10'^ 

Kanupp 

7.24 

X 

10"^ 

-6.0857 X 10'^ 

Pickering 1-8 

8.41 

X 

10-^ 

-5.3829x10'^ 

Point Lepreau 1 

8.22 

X 

lO'^ 

-4.1739x10'^ 

Wolsong 1-2 

8.15 

X 

lO"* 

-4.0335 X 10‘^ 

Atucha 1 

1.25 

X 

10'^ 

-5.4412x10'^ 

Atucha 2 

1.20 

X 

10'^ 

-6.0289x10'^ 
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4.3 Simultaneous Variation of Parameters 


As we have seen above, in each of the reactors studied, bifurcations occur for possibly unrealistic 
values if the parameters other than a,, are considered constants. In this section we vary other 
dimensionless parameters in certain ranges, with a view to find those sets of parameters for 
which bifurcations may occur for values of a,; more likely to occur in practice. The 
dimensionless parameters due to xenon are not considered for the random parameter variation as 
it was seen from Table 13 that the effect of xenon on dynamic bifurcations is very small and 
hence the effect of xenon is not considered in future analysis. For this purpose the values attained 
by the dimensionless parameters for the 220 MWe reactor are treated as reference values and 
denoted by subscript R. These are reproduced in Table 15 along with the ranges explored for 
each parameter. The ranges are selected with the help of Table 6, with a view to allow for the 
worst-case combination. The actual values tried as given in Table 15, are broader than those 
suggested by Table 6. 
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Table 15 ; Ranges of parameter ratios allowed for the worst case PHWR 


Reference Value' 


Range 

Parameter Ratio 

Lower limit 

Upper limit 

3® 

-2.64 

Zfl afR 

1.0 

0.25 


-2.65 

^ ^ ^mR 

2.0 

0.2 


0.0075 

b/bR 

0.3 

3.0 

Cr 

1.08606 

C /Cr 

0.5 

1.25 

9r 

0.02779 

G/Or 

0.5 

1.25 

Pr 

0.00580 

P/PR 

0.01 

10 

q/p 

0.00119 

q/p 

0.2 

5.0 

r /p 

0.00500 

r/p 

0.2 

5.0 


t Reference values are for 220 MWe as given in Tables 2 and 6 and are here denoted by a 
subscript R . 


The worst-case combination of parameters is found using random trials. The results are shown in 
Tables 16 and 17 for static and Hopf bifurcations respectively, together with the corresponding 

values of a*c ■ 
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Table 16: Values of the parameter ratios for the worst case PHWR for the static bifurcation 
within the ranges/ limits given in Table 15 


Parameter Ratio 

Worst-case value for the static bifurcation 


0.25 


-0.2 

b/ba 

3.0 

c/cr 

1.25 

e/ea 

0.5 

P/PR 

10.0 

q/p 

5.0 

r /p 

5.0 

a*c (Core) 

1.7989 X 10*^ K'^ 

a*(. (Core + Moderator) 

1.7534 X 10"^ K'* 

a*(; (Core + Moderator + SG) 

2.5xlO'^K'^ 
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Table 17: Values of the parameter ratios for the worst case PHWR for the dynamic bifurcation 
within the ranges / limits given in Table 15 


Paxameter Ratios 

Worst-case value for the dynamic bifurcation 


P 

3 

11 

O 

a„ = +ve 

am = -ve 

af/a^R 

0.25 

0.25 

0.25 


— 

- 0.2 

2.0 

b/ba 

3.0 

3.0 

3.0 

c/Cr 

0.570 

0.5 

1.333 

G/Sr 

1.25 

0.5 

0.5 

p/PR 

0.01 

0.01 

10.0 

q/p 

— 

5.0 

0.2 

a*c ( K'^) (one-group) 

- 4 . 107 x 10 *^ 

- 6.068 X 10 '^ 

+ 2.122 X 10 "* 


In the worst-case the values of the coolant temperature coefficient of reactivity for which a static 
bifurcation is more likely to occur are found to be of the order of +10-^ K'". While this value is 
physically realistic for PHWRs, It must be noted that this applies when the parameters of all 
reactors are combined in the worst manner. Table 17 also shows the effect of moderator feedback 
on dynamic bifurcation. Even though the physical value of moderator temperature coefficient of 
reactivity (a J for all PHWRs in the present work is taken to be negative, we have allowed a 
small positive value for the moderator temperature coefficient while exploring the ranges for the 
dimensionless parameters, as can be seen horn Table 15. From Table 17 it is clear that the 
worst-case parameters depend strongly on the moderator temperature coefficient of reactivity. 
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The worst-case values of in case of dynamic bifurcation range from -4 x 10'^ to +2 xlO^ K■^ 
depending on the moderator temperatue coefficient of reactivity. 


4.4 Trajectories, Bifurcation Diagrams and 
Phase Portraits 


In this section , we present numerically simulated trajectories and bifurcation diagrams for the 
worst case and for the 220 and 500 MWe. Furthermore, as has been mentioned in Section 3.4 the 
static bifurcation is possible only when the higher-order feedback effects are also present. In the 
figures that follow, these effects are included for the transcritical bifurcation shown. In the case 
of Hopf bifurcation, we have ignored the higher-order feedback effects. 

Figures 5 to 7 shows the time- response for the dimensionless power, and fuel and coolant 
temperatures subsequent to the transcritical bifurcation. The parameter values for which the 
figures are plotted are shown on the figures. It can be seen that while the behavior of the linear 
system, or the system with linear feedback, is divergent as expected, the behavior becomes 
bounded, as soon as any of the higher-order feedback effect is included. The presence of the new 
equilibrium point is indicated by the asymptotic approach to the same m Figures 5 to 7. 

The actual values of the new equilibrium point depend on how far the parameter is from the 
bifurcation point. This is shown in Figures 8 to 10, where 


Epsilon s 



Tlie of stabUity betwoen the operating point and the new fixed point at the biaication 

is also clearly visible in the figures. The soUd curves represent stability while the dashed curves 

represent instability. 
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It should however be mentioned that new equilibrium point is only locally stable ( in accordance 
with the local bifurcation theory). Thus all disturbed situations from the original operating point 
need not be in the basin of attraction of the new fixed point. The behavior shown in Figures 5 to 
7 is valid only for those disturbances which are in the basin of attraction of the new equilibrium 
point. For other disturbances, the reactor behavior will re main divergent. 

Figures 11 to 14 show the behavior subsequent to Hopf bifurcation in the worst-case PHWR 
dynamical system using one group of delayed neutrons. As has already been mentioned , the 
values of the parameters for which this behavior occurs appear to be impractical. These figures , 
however confirm the prediction of the theory predicted in Chapter 3. Figures 15 to 18 show the 
behavior subsequent to Hopf bifurcation in 220 and 500 MWe PHWRs for the core, core along 
with moderator, and core along with moderator and SG, using one group of delayed neutrons. In 
figures only two curves are seen as the curve for the core, and core with moderator overlap. We 
reemphasize that for 220 and 500 MWe reactors, the values of a*<- (Table 14) for which the Hopf 
bifurcations are even more unphysical than those for the worst case (Table 17). However, even in 
this unlikely situation, while the dimensionless power variation during the limit cycles oscillation 
is large, its effect on fuel and coolant temperature is quite limited (1 to 5% for 200 and 500 MWe 
reactors, and ~ 1 5% when the parameters are combined in the worst manner for all the reactors 
listed in Table 6). 

The stability of the limit cycles represented in Figures 11 to 14 is further confirmed by the 
calculation of stability parameter, a, and the characteristic multipliers shown in Table 18. 

We mentioned here that as a > 0 , the situation similar to that shown in Figure 3c occurs in these 
models. In view of the unphysical nature of a*c , we have not further investigated the 
coexistence of stable and unstable limit cycle in the vicinity of bifurcation point. This and its 
associated effects are reported in detail by Manmohan (1996) in the contexts of PWRs and 
BWRs. 
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Table 18; The stability parameter o(0)and the characteristic ( Floquet ) multipliers for the 
PHWR core limit cycle (worst case) 


Parameter 

One Group 

Six Groups 

Stability Parmeter,a(0) 

3.191 xlO"^ 

1.632x10'^ 

Floquet Multipliers 

(e = 0.2) 

0.909 

1.09 X lO'^ 

3.30 X 10"^ 

0.832 

7.14 X 10'^ 
1.76x10'^ 

1.28 X 10'’ 

9.33 X 10'® 
3.48 X 10'® 
1.73x10'® 
1.13x10'® 


4.5 Critical Values for the Control Reactivity 


As given in Section 2.10 , the control may be based on the fractional change in reactor power and 
/ or a suitably defined fractional change in reactor coolant temperature. Table 19 gives the critical 
values, a , , for different values of Kp and tp , as compared to the situation without control (Kp = 
0 , tp = 0). It can be seen that bifurcation becomes all the more unlikely for all values of Kp < 0, 

the only situation that can occur in practice. 
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Table 19. Critical values of the coolant temperature coefficient of reactivity (aj for different 
values of Kp and tp for the Reference Reactor (220 MWe ; Core + moderator + SG) using 6 
groups of delayed of neutrons 




Critical values, , K'^ 

Kp 


Static 

Dynamic 

0 

— 

2 X 10'^ 

-6.11x10'^ 

-0.1 

0 

2x10'^ 

-1.71 X 10^ 


0.1 

3.53 xlO"* 

-1.078 

-10 

0 

2 X 10'^ 

- 1.53 X 10^ 


0.1 

3.35x10'^ 

-1.02x10^ 


Table 20 shows the critical values of for different values of t;. . It is seen that K is of the 
order of -1. This, for the values of T^q, Tjo and p given in Tables 1 and 3 amounts to a coolant 
reativity feedback of the order of -10 $ per degree rise of D 2 O coolant temperature. A practical 
value for the same which is likely to be of the order of -10 ^ to 10 $ per degree rise in coolant 
temperature. Thus it can be concluded that control system is unlikely to introduce any kind of 
bifurcation. 
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Table 20 : Critical values of the control parameter K, ( Kp = 0 , tp = 0 ) for Hopf bifurcation in 
Reference PHWR(220 MWe) 



♦ 

tc 

Kc 

0 

-1.22 

0.001 

-1.24 

0.01 

-1.40 


4.6 Conclusions And Observations 


As summarized in the preceding sections, an analysis of Hopf bifurcation and static (transcritical) 
bifurcation in some models of PHWR dynamical systems has been presented in this work. Based 
on the analytical and numerical results , the following conclusions are drawn: 

• Hopf and transcritical bifurcations in the models considered occur in one-dimensional 
parameter space, i.e., it is of codimension one. In general, the critical values of the variable 
parameter calculated using different mathematical models are found to differ considerably 
from each other. This is particularly so if the models based on the effective-lifetime and 
prompt-jump approximation are used and hence these models should not to be used for 
bifurcation studies. It is pointed out in this work that qualitatively, the behavior of the models 
with six groups of delayed neutrons is similar to that of the one-group models. 

• In lumped parameter PHWR dynamical systems based on the two temperature feedback 
model and the parameter ranges explored, it appears that a PHWR core may be an unlikely 
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candidate to experience a Hopf bifurcation instability. This remains true for all (negative) 
gains in the proportional controllers and any time constants associated with the differential 
controllers. Effect of xenon is seen to be very low and hence its effect is omitted in the 
analysis of the present work. However, when the moderator feedback is included, for 
negative moderator temperature coefficient of reactivity am worst-case analysis shows that 
the likelihood of Hopf bifurcation is somewhat increased. This is not the case for the original 
three-temperature feedback PHWR model for any of the reactors. 

• Static or transcritical bifurcation can occur only if the higher-order feedback coefficients 
are present. Of the three types of static bifurcation , it is shown in this work that oni> 
transcritical bifurcation can occur. The static bifurcation occurs for the positive values of a.^ 
and is of the order of 10'^ to lO'^ K'‘, while the actual values is of the order 10'^ to 10 K . 
For the worst-case, the critical value of for the core along with the moderator and SG is of 
the order of 10'^ K"'. However this is so if all the parameters are combined m the worst 
possible manner, an unlikely situation to occur in practice. 

• It is noted in the present work that subsequent to transcritical bifurcation, the new 
equilibrium point is only locally stable and may not be an attractor for the all disturbances 
from the normal operating condition. Thus there remain disturbances of the operatmg 
condition from which the reactor will show unstable divergent behavior mspite of the 
presence of higher-order feedback coefficients. 

A few observations and suggestions which may be useful in future extensions of the pres 

I 

work are given below; 

. We have used space-independem lumped parametet models for the present study. It would be 
desirable to extend this study to include spafral effects. 
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• In present work, the control reactivity whether based on the change in reactor power or based 
on the change in coolant temperature, is assumed to act instantly. In future studies, 
appropriate governing equations for the control reactivity insertion should be included. 

• The present study has been restricted to long-term or asymptotic behavior in the dynamical 
systems studied. The role of nonlinearities in the transients, which may sometimes last veiy 
long and whose behavior may be important in fission reactor dynamics and safety, should be 
include in future studies. This may also have an important bearing on the optimal design of 
the controllers / control law. 

• In the present work, we have studied the local behavior in the vicinity of the bifurcation 
points. A global analysis may be attempted in future studies. 

• Nonlinearities also effect the transient behavior of the reactor under normal operating 
parameter values. Their effect on the transients may have an important bearing on the optimal 
design of the controllers / control law. This can be a fruitful exercise to undertake in future. 
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Figure 5: Time variation of dimensionless power p subsequent to transcntical bifurcation m a 
PHWR core (worst case) 
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Figure 6 : Effect of Steam Generator on time variation of dimensionless fuel 

temperature Oy subsequent to transcritical bifurcation using one-group of delayed neutrons in 

PHWR core with moderator (worst case) 

s=i5,^=^^omK-^ 

[ aj- ac J 
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Figure 7; Effect of one and six groups of delayed neutrons on time variation of dimensionless 
fuel temperature 3y subsequent to transcritical temperature in PHWR core with moderator 



i.o,^=^ = o.oir"’ 
«/ CLc J 
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Figure 1 1 : Behavior of Dimensionless Power (neutron flux ) versus time in PHWR core 
(worst case) subsequent to Hopf bifurcation using one-group of delayed neutrons {s = 02) 
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Figure 12 : Behavior of Dimensionless coolant temperature (neutron flux ) versus time in PHWR 
core (worst case) subsequent to Hopf bifurcation using one-group of delayed neutrons = 02) 
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Figure 13 : Phase portrait of limit cycles in PHWR core with moderator(worst case) on p - 5 
plane subsequent to Hopf bifurcation = 0.2) 
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Figure 14 ; Phase portrait of limit cycles in PHWR core with moderator(worst case) on 3 
plane subsequent to Hopf bifurcation ( ^ = 0.2 ) 
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Figure 15 : Phase portrait of limit cycles in PHWR (220 MWe) for the core, core along with 
moderator, and core with moderator and Steam Generator on - 3^ plane subsequent to Hopf 
bifurcation using one-group of delayed neutrons(e = 0.2 ) 
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Dimensionless Power 


Figure 16 ; Phase portrait of limit cycles in PHWR (220 MWe) for the core, core along with 
moderator, and core with moderator and Steam Generator on p - 3^ plane subsequent to Hopf 
bifurcation using one-group of delayed neutrons(£ = 0.2) 
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Figure 17 ; Phase portrait of limit cycles in PHWR (500 MWe) for the core, core along with 
moderator, and core with moderator and Steam Generator on - 3^ plane subsequent to Hopf 

bifurcation using one-group of delayed neutrons(5 = 0.2 ) 
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Figure 18 : Phase portrait of limit cycles in PHWR (500 MWe) for the core, core along with 
moderator, and core with moderator and Steam Generator on p — plane subsequent to Hopf 
bifurcation using one-group of delayed neutrons(ff = 02) 
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Appendix A 


Sample calculations of Cf , Co and 

Estimation of Cf 


Cf = Cf nif = 17.92 MJK'\ where 

specific heat of fuel elements, Cf = 0.32 KJkgK'^ 

fuel inventory, mf = 56 x 10^ kg 

Vf = mf / Pf = 56000 / 11000 = 5.09 m^, where Vf is volume of fuel elements, and 
Pf= 1 1000 kgm' is the density of fuel material. 

Estimation of Cp 

fuel pellet diameter , d = 1.427 cm 
clad thickness , t = 0.041 cm 
Vd«/V,„, = 4t(d + ty <1^ = 0.118 

V (total volume of 306 pressure tubes/ coolant channels) = (ti / 4) d; Lx 306 

= (71 / 4) (0.08255)^ X 6 X 306 
= 9.89 m^ 
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where 

inner diameter of pressure tube , d-, = 0.08255 .m 
active length of the pressure tube , L = 6 m ' • 

Vf/V = 5.09/ 9.83 = 0.518 

Vciad / V = / Vf„ei) X (Vr,„/ V) = 0.0613 

Vsm.cture / V = / V ) X 1 .2 = 0.0735 

(20% additional structure volume over and above the clad volume allowed in the pressure tubes.) 

(y fuel ^Structure )/Vt„,^ = 0.5915 
Voooi«nt/V = 0.4085 s 0.41 
■Vcooiant = 0-41 X 9.83 = 4.015 m^ 

^coolant ^coolant X pc " 3.21' toones, where 
density of coolant , = 800 Kg m' 

Tliis value of mediant has ben approximated to 3 tonnes. 

Cc = iBe Cg = 16.5 MJK'*, where 
specific heat of coolant , c^ = 5.5 KJkg K‘^ 

EsVsjTiation of 

Cm = Ittm Cm. wlieie 

is mass of moderator in the calandxia, and 
specific heat of the moderator , c^, = 4.2 kJK 

has been estimated as follows: 

Calandria volume = ( it / 4) L s 100 m^. 


8S 


where D and L axe the diameter and active length of the calandria. 

For the present calculation D = 5m , L s 5m. 

Volume displaced by the coolant channels / pressure tubes = ( tt / 4) d^^ L x 306 

= (7t/ 4) (0.12)^x5x306 
= 17.3 m^ 

where 

outer diameter of the pressure tube , d^ s 0.0012 cm 

Vmoderator = (100.0 - 17.3 ) m^ = 82.7 m' 
mm = Vmoderator X Pm = 80 tonnes, where 
density of the moderator, = 1 100 kg m'^ 

Thus approximated 60 % of the total moderator calculation inventory of 137 tonnes (NPCIL, 
1992) is in the calandria. We have made a similar assumption for the other reactors to calculate 
the mass of the moderator in the calandria. 


89 



